%% load data
function [next_trace_cell] = ...
    load_fun(pstruct,...
        TT_datetime_local,...
        main_starttime,...
        batchin)
    %% Download on one core at time to avoid problems
    %josh crozier 2020
    next_trace_cell = cell(pstruct.events_per_file+1,1);
    for eventin = 1:pstruct.events_per_file
        try
%             if pstruct.useirisfetch
%                 %% download data with irisfetch (not working currently)
%                 %args: Net, Sta, Loc, Cha, Starttime, Endtime [,quality][,includePZ][,verbosity]
%                 % D � The state of quality control of the data is indeterminate.
%                 % R � Raw Waveform Data with no Quality Control
%                 % Q � Quality Controlled Data, some processes have been applied to the data.
%                 % M � Data center modified, time-series values have not been changed.
%                 traces=irisFetch.Traces(pstruct.networks, pstruct.stations,...
%                     pstruct.locations, pstruct.channels,...
%                     datestr(TT_datetime_local(eventin) - pstruct.startshift, 31),...
%                     datestr(TT_datetime_local(eventin) + pstruct.endshift, 31),...
%                     'includePZ');
% 
%                 elapsedtime = datetime - main_starttime;
%                 disp(strcat('batch:',num2str(batchin),' event:',num2str(eventin),...
%                     ' tracelen:',num2str(length(traces)),' runtime=',string(elapsedtime)));
% 
%                 if ~isempty(traces)
%                     %combine traces
%                     traces = combine_traces(traces);
% 
%                     %check for traces not spanning full duration
%                     late = datetime([traces.startTime],'ConvertFrom','datenum')...
%                         - (TT_datetime_local(eventin) - pstruct.startshift);
%                     early = datetime([traces.startTime],'ConvertFrom','datenum') ...
%                         + seconds([traces.sampleCount]./[traces.sampleRate]) ...
%                         - (TT_datetime_local(eventin) + pstruct.endshift);
%                     incompind = late > 2*seconds([traces.sampleRate].^-1) | ...
%                         early < -2*seconds([traces.sampleRate].^-1);
%                     if any(incompind)
%                         strcat('incomplete traces:', ...
%                             datestr(TT_datetime_local(eventin)-pstruct.startshift),31)
%                         traces = traces(~incompind);
%                     end
%                 end
% 
%                 if ~isempty(traces) 
%                     %check for bad data
%                     tracebad = false(size(traces));
%                     for tracein = 1:length(traces)
%                         if any(isnan(traces(tracein).data))
%                             tracebad(tracein) = true;
%                         end
%                     end
%                     if any(tracebad)
%                         strcat('bad data:', ...
%                             datestr(TT_datetime_local(eventin)-pstruct.startshift),31)
%                         traces = traces(~tracebad);
%                     end
%                 end
% 
%                 if ~isempty(traces)
% 
%                     %remove responses
%                     traces = remove_response(traces, pstruct.water_level);
%                 end
%             else

            %% GISMO download or load
            starttemp = TT_datetime_local(eventin) - pstruct.startshift;
            endtemp = TT_datetime_local(eventin) + pstruct.endshift;
            if ~pstruct.uselocal
                %iris download
                wf = waveform(pstruct.ds, pstruct.ctags, ...
                    datestr(starttemp, 31),...
                    datestr(endtemp, 31));
                
                %attach local responses (due to iris bug)
                for wfin = 1:length(wf)
                    if ~isempty(wf(wfin))
                        temp_channel_station_str = strcat(get(get(wf(wfin),'ChannelTag'),'network'),...
                            "_",get(get(wf(wfin),'ChannelTag'),'station'),...
                            "_",get(get(wf(wfin),'ChannelTag'),'channel'));
                        chanin = find(strcmp(temp_channel_station_str,pstruct.channel_station_str),1);
                        temp_sacpz = [];
                        for respin = 1:length(pstruct.resps{chanin})
                            %check if waveform starts before end of each epoch
                            if isempty(pstruct.resps{chanin}(respin).endtime)
                                beforeend = true;
                            elseif datenum(get(wf(wfin),'Start')) < pstruct.resps{chanin}(respin).endtime
                                beforeend = true;
                            else
                                beforeend = false;
                            end
                            %check if waveform also starts after beginning of each epoch
                            if beforeend &&...
                                    datenum(get(wf(wfin),'Start')) > pstruct.resps{chanin}(respin).starttime
                                temp_sacpz = pstruct.resps{chanin}(respin);
                                break
                            end
                        end
                        if ~isempty(temp_sacpz)
                            wf(wfin) = addfield(wf(wfin),'SACPZ',temp_sacpz);
                            set(wf(wfin),'SACPZ',temp_sacpz);
                        end
                    end
                end
                
            else
                %load from local sac files
                wf = gismo_load_sac(pstruct.localdrive,...
                    starttemp,endtemp,...
                    pstruct.ctags,pstruct.responses);
            end
            
            for win = 1:length(wf)
                if ~isempty(wf(win))
                    %pad in case of gaps at ends
                    wf(win) = pad(wf(win),datenum(starttemp),datenum(endtemp),NaN);
                end
            end

            elapsedtime = datetime - main_starttime;
            disp(strcat('batch:',num2str(batchin),' event:',num2str(eventin),...
                ' tracelen:',num2str(length(wf)),' runtime=',string(elapsedtime),...
                ' start:',datestr(TT_datetime_local(eventin) - pstruct.startshift,31)));

            %check for bad waveforms
            if ~isempty(wf)
                keepind = true(size(wf));
                for win=1:length(wf)
                    if isempty(wf(win))
                        keepind(win) = false;
                    else
                        try
                            get(wf(win),'SACPZ'); %some are missing response data
                        catch
                            keepind(win) = false;
                            disp(strcat('noresp_',get(wf(win),'station'),'_N=',num2str(eventin)))
                        end 

                        %toss flat waveforms
                        if range(get(wf(win),'data')) == 0
                            keepind(win) = false;
                        end

                        %toss waveforms with big gaps
                        if keepind(win)
                            gapind = ~isfinite(get(wf(win),'data'));
                            gaplen = floor(seconds(pstruct.gapdur)*get(wf(win),'freq'));
                            if sum(gapind) >= gaplen
                                gapcount = gapind;
                                for ii = 1:gaplen
                                    gapcount = gapcount + [zeros(ii,1);gapind(1:end-ii)];
                                end
                                if any(gapcount >= gaplen)
                                    keepind(win) = false;
                                end
                            end
                            clear gapind
                        end
                    end
                end %end loop over waveforms
                wf = wf(keepind);
            end

            if ~isempty(wf)
                %fill remaining gaps
                wf = fillgaps(wf,'interp');

%             % Deconvolve instrument response
%             wf = detrend(wf);
%             filterObj = filterobject('b',[smin/4, 4*smax],2);
%             wfFilt = filtfilt(filterObj,wf);
%             wfCorrected = response_apply(wfFilt,filterObj,'antelope',antelopeds);

                % convert to traces struct
                traces = struct('network', get(get(wf,'ChannelTag'),'network'),...
                    'station', get(get(wf,'ChannelTag'),'station'),...
                    'location', get(get(wf,'ChannelTag'),'location'),...
                    'channel', get(get(wf,'ChannelTag'),'channel'),...
                    'data', get(wf,'data'),...
                    'sampleRate', num2cell(get(wf,'freq')),...
                    'startTime', num2cell(get(wf,'start')),...
                    'sacpz', get(wf,'SACPZ'));%,...
%                     'latitude', num2cell(get(wf,'LATITUDE')),...
%                     'longitude', num2cell(get(wf,'LONGITUDE')),...
%                     'elevation', num2cell(get(wf,'ELEVATION')));%,...
%                     'sensitivity', num2cell(get(wf,'SENSITIVITY')),...
%                     'sensitivityFrequency', num2cell(get(wf,'SENSITIVITYFREQUENCY')),...
%                     'sensitivityUnits', num2cell(get(wf,'units')),...
%                     'instrument', num2cell(get(wf,'INSTRUMENT')),...
%                     'depth', num2cell(get(wf,'DEPTH')),...
%                     'azimuth', num2cell(get(wf,'AZIMUTH')),...
%                     'dip', num2cell(get(wf,'DIP')));
                endTime = num2cell([traces.startTime] + [get(wf,'duration')]);
                [traces.endTime] = endTime{:};
                for tracein = 1:length(traces)
                    traces(tracein).sampleCount = length(traces(tracein).data);
                    traces(tracein).startTime = ...
                        datetime(traces(tracein).startTime,'ConvertFrom','datenum');
                    traces(tracein).endTime = ...
                        datetime(traces(tracein).endTime,'ConvertFrom','datenum');
                    if ~pstruct.uselocal
                        traces(tracein).latitude = (get(wf(tracein),'LATITUDE'));
                        traces(tracein).longitude = (get(wf(tracein),'LONGITUDE'));
                        traces(tracein).elevation = (get(wf(tracein),'ELEVATION'));
                    end
                end
                if pstruct.uselocal & pstruct.use_synthetics
                    [traces,keepind] = local_HV_locations(traces);
                    traces = traces(keepind);
                    %wf = wf(keepind);
                end
                clear wf
                
%                 if range([traces.startTime]) > seconds(1) || range([traces.endTime]) > seconds(1)
%                     'debug pause'
%                 end
                
                if pstruct.use_synthetics
                    traces = make_mogi_synthetic_waveform_fun(traces,...
                        pstruct.synthetic_T,...
                        pstruct.synthetic_Q,...
                        pstruct.synthetic_A,...
                        pstruct.synthetic_step,...
                        pstruct.synthetic_noise,...
                        pstruct.synthetic_sourcenoise,...
                        pstruct.synthetic_onset+pstruct.startbounds(1),...
                        pstruct.synthetic_mu,...
                        pstruct.synthetic_nu,...
                        pstruct.synthetic_lat,...
                        pstruct.synthetic_lon,...
                        pstruct.synthetic_elevation);
                end

                %remove responses
                traces = remove_response(traces, pstruct.water_level);
                
                %smooth + downsample
                Ts = pstruct.Tmin/6;
                span = Ts;
                time = traces(1).startTime + ...
                    (seconds(0):Ts:seconds(1./traces(1).sampleRate*traces(1).sampleCount))';
                traces = interp_traces(traces, time, span, pstruct.taperfrac);
                
            else
                traces = [];
            end

            if ~isempty(traces)
                next_trace_cell{eventin} = traces;
            else
                next_trace_cell{eventin} = [];
            end
            clear traces

        catch ME
            %check for downloading errors
            disp(getReport(ME))
            disp(strcat('error downloading:',num2str(batchin),',',...
                num2str(eventin),',',datestr(TT_datetime_local(eventin),31)))
            next_trace_cell{eventin} = [];
%                     'rethrowing errors'
%                     rethrow(ME)
        end

    end %end downloading loop

end %end download function


